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Abstract. The development of parallel supercomputers allows today the detailed study of the collapse and the frag- 
mentation of prestellar cores with increasingly accurate numerical simulations. Thanks to the advances in sub-millimeter 
observations, a wide range of observed initial conditions enable us to study the different modes of low-mass star forma- 
tion. The challenge for the simulations is to reproduce the observational results. Two main numerical methods, namely 
AMR and SPH, are widely used to simulate the collapse and the fragmentation of prestellar cores. We compare here thor- 
oughly these two methods with numerical resolution requirements deduced from previous studies. Our physical model is 
as simple as possible, and consists of an isothermal sphere rotating around the z-axis. We first study the conservation of 
angular momentum as a function of the resolution. Then, we explore a wide range of simulation parameters to study the 
fragmentation of prestellar cores. There seems to be a convergence between the two methods, provided resolution in each 
case is sufficient. Resolution criteria adapted to our physical cases, in terms of resolution per Jeans mass, for an accurate 
description of the formation of protostellar cores are deduced from the present study. This convergence is encouraging 
for future work in simulations of low-mass star formation, providing the aforementioned criteria are fulfilled. 

1 Introduction 

Star formation is known for being the place of extreme variations in length and density scales. Although it is established 
that stars form in dense cores, the non-linear evolution makes it difficult to perform accurate calculations of the collapse 
and the fragmentation of a prestellar core. The star formation process is the outcome of complex gas dynamics involving 
non-linear interactions of gravity, turbulence, magnetic field and radiation. Klein et al. (2007) point out that developing 
a theory for low-mass star formation remains one of the most elusive and important goals of theoretical astrophysics. 
Numerical simulations allow today to study star formation processes and, thanks to the recent development of super 
computers, a significant increase in the description of the dynamical range of low-mass star formation has been reached. 
The computational challenge stems from the fact that star formation occurs in clouds over many orders of magnitude 
in spatial and density scales. The main issue for the simulations is to follow the gravitational collapse while resolving 
precisely the Jeans length, which scales as Aj oc p~ 1 / 2 for an isothermal gas. 

Different approaches are used to study star formation through numerical simulations and include more and more 
detailed physics. One key question resides in the choice of the most appropriate numerical method to be used to study low- 
mass star formation. Nowadays, two completely different approaches can handle this problem with sufficient accuracy: the 
Adaptive Mesh Refinement method for Eulerian grids (AMR) and the Smoothed Particle Hydrodynamics method (SPH) 
for a Lagrangian approach. We use the AMR code RAMSES (Teyssier 2002) and the SPH code DRAGON (Goodwin et 
al. 2004). The debate about whether one method is most appropriate remains open and no systematic comparison has been 
done in the domain of low-mass star formation. Recently, in a comparative study between SPH and AMR, Agertz et al. 
(2006) have shown that the SPH method is of limited accuracy for describing Kelvin-Helmotz instabilities in the presence 
of a strong density gradient. This effect in SPH simulations is due mainly to the cyclic, kick-drift phenomenon, which 
strongly depends on the form of the SPH force calculation coupled with the treatment of neighbor particles. In the low- 
mass star formation field, Fromang et al. (2006) compares quite successfully AMR hydrodynamical collapse calculations 
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with SPH ones. Several studies provide key starting points for the numerical investigation of low-mass star formation. 
Truelove et al. (1997) give an empirical criterion for the Jeans length resolution in AMR calculations to avoid spurious 
numerical fragmentation while Bate & Burkert (1997) provide a minimum resolution criterion for SPH calculations with 
self-gravity to correctly model fragmentation. In the present study, we compare thoroughly the two approaches in the 
context of low-mass prestellar core formation. We focus on the dependency of angular momentum conservation and 
fragmentation on physical and numerical initial conditions, in order to derive resolution criteria adapted to our study. 



2 Definitions of the test cases 



To make comparison between codes easier, we adopt simple initial conditions, similar to those chosen in previous studies 
(Bate & Burkert 1997). We consider an uniform sphere of molecular gas of initial radius Rq, rotating around the z-axis 
with a constant angular velocity Q . We set up the cloud mass at Mq = 1 M Q and the temperature at 10 K. For a 
mixture of molecular hydrogen, helium and heavy elements, with mean molecular weight p = 2.2, this corresponds to an 
isothermal sound speed of Co ~ 0.19 km.s~'. For the case where fragmentation occurs, we use a m = 2 azimuthal density 
perturbation. The initial energy balance of our model is determined by two dimensionless parameters corresponding to 
the ratio a between the thermal energy and the gravitational energy and to the ratio f3 of the rotational and the gravitational 
energy: 

5 R kT 1 #o Q o 

a= :B= — — . (2.1) 
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Since we use a constant initial mass of 1 M and a constant temperature, changing one of the two parameters, namely a, 
gives the sphere radius Rq. The higher a, the larger Rq. The angular velocity is given by the parameter f3. 

In order to mimic the thermal behaviour of a star-forming gas, we use a barotropic equation of state. Masunaga & 
Inutsuka (2000) showed that the core follows closely a barotropic equation of state, providing a good approximation 
without resolving radiative transfer. We use 
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where Co is the isothermal sound speed at 10 K and p c = 10~ 13 g.cm~ 3 is the critical density which corresponds to the 
transition from an isothermal to an adiabatic state (Larson 1969). At low densities, p «c p c , C s ~ Co = 0.19 km.s~'. The 
molecular gas is able to radiate freely by coupling thermally to the dust and therefore remains isothermal at 10 K. At high 
densities p > p c , we assume that the cooling due to radiative transfer is trapped by the dust opacity. Therefore, P oc p 5 ^ 3 
which corresponds to an adiabatic monoatomic gas with adiabatic exponent y = 5/3. Note that molecular hydrogen 
behaves like a monoatomic gas until the temperature reaches several hundred Kelvin, since the rotational degrees of 
freedom are not excited at lower temperatures, and hence y = 5/3 is the appropriate adiabatic exponent (e.g. Masunaga 
& Inutsuka 2000). 



3 Numerical methods and refinement criteria 



3. 1 AMR code RAMSES 

In this paper, we use the AMR code RAMSES (Teyssier 2002), which integrates the "tree-based" data structure allow- 
ing recursive grid refinements on a cell-by-cell basis. RAMSES combines a tree-based AMR grid and a second order 
Godunov hydrodynamical scheme coupled with a gravity solver. The Godunov hydrodynamical solver is able to capture 
discontinuities with a high precision level. The equations solved in RAMSES are the Euler equations in their conservative 
form. Furthermore, it has the possibility to use variable timesteps at each refinement level. Concerning time integration, 
RAMSES uses a second-order midpoint scheme, where positions and velocities are updated by a predictor-corrector step. 
Recently, an ideal MHD version of RAMSES has been developed by Fromang et al. (2006). 

Our refinement criterion is based on the Jeans length resolution which is necessary to treat accurately gravitational 
collapse. We impose a minimum number of points Nj per Jeans length A 3 . The cells dimension must be smaller than a 
constant fraction of the local Jeans length. The dimension of cells belonging to the ( t refinement level is Lb ox /2 (i , where 
Lbox is the physical length of the simulation box. The mesh is locally refined in order to verify the local Jeans criterion: 
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Truelove et al. (1997) defined a minimum resolution condition for the validity of grid-based simulations aimed at 
modeling the collapse of a molecular cloud core, namely A^j > 4. This condition ensures that the collapse is of physical 
rather than of numerical origin. 



3.2 SPH code DRAGON 

We use the SPH code DRAGON (Goodwin et al. 2004). The SPH formalism relies on an interpolation method which 
allows any function to be expressed in terms of its values at the location of various points called particles. For numerical 
works, the integral interpolant for the variable A(r,) is approximated by a summation interpolation over the particle's 
nearest neighbors: 



t i 
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where Aj is the value associated with particle j, h/j = (hi + hj)/2 and hi is the adaptive smoothing length of particle 
i, defined such that the particle kernel volume (i.e the resolution element) contains a constant mass, i.e. a constant 
number of neighbors N^. The main advantage of SPH is its Lagrangian conservation property in contrast with grid- 
based methods. Resolution elements are then concentrated in high density regions in SPH methods whereas the AMR 
allows high resolution of all regions in the flow. The standard SPH formalism adopts artificial viscosity. Some alternative 
formalism such as Godunov SPH has been proposed in order to avoid the use of artificial viscosity, but these methods are 
not yet mature. 

In the SPH, the resolution in mass is fixed and thus the Jeans length resolution deteriorates with increasing density 
during the isothermal phase, contrary to the AMR. Bate & Burkert (1997) showed that the behaviour of a Jeans-mass 
clump of gas with radius ~ h is dominated by the numerical implementation. The minimum resolvable mass must then be 
larger than the interpolation mass. Bate & Burkert (1997) take the smallest mass than can be resolved in SPH calculations 
to be equal to the mass of ~ 2A^j particles. We can determine our initial number of SPH particles according to this 
criterion. The Jeans mass is Mj ~ 6G~ 3/ V 1/2 C S 3 and the minimum resolvable mass is M res = niN^, where m is the 
particle mass. Hence, we can define a Jeans condition corresponding to the minimum value of C^p -1 ^ 2 , given by the 
barotropic equation of state d2.2l ). i.e. 2 3 / 2 C 3 ) Pr '^ 2 : 

2 3/2 6C 3 5.35X10- 3 
m < m max — M . (3.3) 

Considering a spherical cloud of mass Mq = 1 M , the initial number of particles N p has to verify N p > Mo/m max ~ 9300 
if A^n = 50. This is the critical number of particles used in SPH calculations to study the collapse of a dense core. We 
have in that case exactly 2N N (i.e. two resolution elements) particles per critical Jeans mass. 



4 Results 

4. 1 Angular momentum conservation 



Table 1. Summary of the different simulations (left table: SPH, right table: AMR) performed to study angular momentum conservation. 

for the AMR calculations gives us the number of cells describing the initial sphere. The quantity N com representing the number of 
cells/particles with density p > 1 x 10~ 15 g.cm" 3 . Time to is the time when the density becomes higher than p c = 1 x 10~ 13 g.cm -3 . 
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We start by comparing the global properties of the collapse in the two codes in the simple case of a uniform sphere 
rotating with a constant angular velocity around the z-axis, with no perturbation. We look at the collapse time, the 
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accretion shock (see Commercon et al. 2007) and finally the angular momentum conservation. In this section, we 
focus on this last issue. Considering our axisymetric model, without azimuthal perturbation, we can easily investigate 
the effect of numerical resolution on local angular momentum conservation, since local angular momentum should be 
perfectly conserved. The loss of local angular momentum in our model is only due to unphysical transport inherent to 
the numerical methods used in the two codes. Thanks to its Lagrangian properties, the SPH calculation gives access to the 
angular momentum that each particle has initially, i.e. the angular momentum that particle should have if the numerical 
scheme was conserving it exactly. Having access to the particle initial angular momentum and to the same quantity at a 
given time, the loss of angular momentum is easily calculated. 

We carried out a set of simulations within a wide range of resolution parameters. The initial sphere has the parameters 
a = 0.65 corresponding to an initial radius 9.2 x 10 16 cm and a density po ~ 6.02 x 10~ 19 g.cirT 3 . The corresponding 
freefall time is tff = (3n/32Gpo)^ 2 ~ 86 kyr. We set f3 = 0.01, corresponding to an orbital time t lot = 2.8 x 10 3 kyr. 
Table Q] summarizes the different SPH and AMR calculation runs for this case. Calculations converged with increasing 
numerical resolution to a value slightly greater than the theoretical freeffall time because of the rotational support. 

Figure Q] (left) displays the azimuthal velocity component as a function of the radius r on the xy-plane for the SPH 
simulations. It is obvious that low resolution simulations are not able to conserve properly the angular momentum. With 
5 x 10 3 particles, we obtain contra-rotating particles at the center. This pure numerical effect shows that a minimum 
resolution is required even for a simple collapse model. It appears that a minimum of 5 x 10 4 particles is required to get 
an acceptable angular momentum conservation, within less than 10%, for the case of the present study. The improvement 
of angular momentum conservation with a larger number of particles eventually saturates for large numbers of particles. 
We checked that using a larger number of neighbours does not improve the conservation of local angular momentum. 




Fig. 1. Azimuthal velocity at to as a function of the radius on the equatorial plane for SPH (left) and AMR (right) calculations at 
corresponding t . The left-hand plot (Fig. QJ) shows SPH results for various N p and a constant = 50. The solid line represents the 
theoretical azimuthal velocity interpolated at t and is denoted as dM/dt = 0. The right-hand plot (Fig. [Tp) shows AMR results for 
Nj = 10 and t^n = 5, 6 and 7. The theoretical azimuthal velocity is plotted also for easy comparison with the SPH results. 

Figure Q] (right) shows results obtained with the RAMSES code. AMR curves are plotted and compared with the 
theoretical one obtained previously for the SPH. The simulations with £ m j n = 6 and £ m \ n = 7 give approximately the same 
results. In all these AMR simulations, the core resolution is better than with the SPH (see parameter A^oie in TableHJ. For 
an initial resolution of Anin=5, the AMR scheme does produce some angular momentum lag. This can be due to the fact 
that with a poor initial resolution, the interpolation of the gravitational potential tends to convert gravitational energy into 
rotational energy. 



4.2 Fragmentation 

This subsection is devoted to the exploration of various numerical parameters in a physical case known to fragment. In 
Commercon et al. (2007), we study the effect of varying the initial grid resolution £ m ; n and the number of cells within a 
Jeans length iVj for AMR calculations. We also present our SPH calculations with various and iVp. We carry out this 
study for three different cases with various thermal support. Here, we focus on the case with initial parameters: a = 0.5, 
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/? = 0.04, giving p = 1.35 X 10~ 18 g.cirT 3 , R = 7.07 x 10 lf> cm, Qo = 2.12 x 10~ 13 rad.s" 1 and t ff = 57 kyr. The initial 
perturbation amplitude is A = 0.1. 

Figure [2] shows density maps in the equatorial plane for the results of two amongst the most resolved calculations at 
three timesteps, namely, from top to bottom, to+5 kyr, to+6 kyr and to+7 kyr. The left column shows maps for AMR 
calculations with C m \ n = 7 and Nj — 15 whereas the right column displays SPH maps for calculations with N p = 5 x 
10 s and Nn = 50. We show here the two converged calculations of each method that we first studied independently. 
The convergence between the two methods for these physical and numerical parameters set is striking for the two first 
timesteps. The calculations give the same fragmentation time and pattern, although satellites and the central object are 
bigger with the SPH. Symmetry is eventually broken with the SPH calculations, whereas the AMR ones preserve the 
symmetry longer. 
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Fig. 2. AMR and SPH calculation density maps in the xy-plane at three different times for the case a = 0.50, /3 = 0.04. The times 
correspond to to + 5 kyr, to + 6 kyr and to+7 kyr, from to bottom, respectively. The AMR calculations plotted on the left column have 
been run with ( mm = 7 and Nj = 15. The right column shows the results for the SPH calculations with N p = 5 x 10 5 and A?n = 50. 



5 Conclusion 



We show that we reach good convergence between AMR and SPH methods provided one uses sufficient numerical re- 
sources. First, we take a simple model to study local angular momentum conservation. The initial study shows that 
angular momentum is better conserved with the AMR approach. As shown in Fig. QJ, a smaller number of particles in 
SPH calculations leads to poor local angular momentum conservation. In AMR calculations, a poor initial computational 
domain resolution (i.e. C m \ n < 6) leads to unphysical transfer of gravitational energy to rotational energy (see Fig. [TJi). 
A significant loss of angular momentum will affect fragmentation since less rotational support can balance gravitational 
collapse. The smallest parameter set for SPH calculations, required to go through gravitational collapse without signifi- 
cant loss of angular momentum, corresponds to a number of ~ 530 particles per Jeans mass at the critical density p c , i.e. 
5 particles per Jeans length. The equivalent minimum resolution criterion for AMR calculations is C m \ n > 6 and Nj = 4. 
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For the fragmentation case, the two approaches show good agreement for the general pictures. Details are better 
resolved in AMR calculations thanks to the refinement method based on the local Jeans length, whereas the resolution 
deteriorates with increasing density with the SPH, because of the fixed mass resolution. Numerical calculations of pro- 
tostellar collapse should thus be conducted with great care, with a detailed examination of numerical resolution. The 
resolution criteria to reach convergence for the cases we explore in details is about ~ 5000 particle per Jeans mass in 
SPH and A^j = 15 in AMR, and could become prohibitive in some cases, particularly for low thermal support. In such a 
case, it is more difficult to get a convergence. The horizon of predictability is very short, since fragmentation occurs very 
quickly. The latter parameter set seems to be a lower resolution limit for dense core collapse and fragmentation with SPH 
calculations. Using a lower number of particles will lead to spurious early fragmentation due to numerical effects. The 
present work can be used to assess the validity of numerical tools to study star formation. 

Calculations have been performed thanks at the PSMN (ENS Lyon) and CCRT (CEA) supercomputating facilities, as well as on the CEMAG computing 
facility supported by the French ministry of research and education through a Chaire d' Excellence awarded to Steven Balbus. 

References 

Agertz, O., Moore, B., Stadel, J., Potter, D., et al 2006, astro-ph/0610051| 
Bate, M. R. & Burkert, A. 1997, MNRAS, 288, 1060 

Commercon, B., Hennebelle, P., Audit, E., Chabrier, C. & Teyssier, R., submitted to A&A 

Fromang, S., Hennebelle, P. & Teyssier, R. 2006, AAP, 457, 371 

Goodwin, S. P., Whitworth, A. P. & Ward-Thompson, D. 2004, AAP, 414, 633 

Klein, R. I., Inutsuka, S.-L, Padoan, P. & Tomisaka, K. 2007, Protostars and Planets V, 99-116 

Larson, R. 1969, MNRAS, 145, 271 

Masunaga, H. & Inutsuka, S -i. 2000, ApJ, 531, 350 

Teyssier, R. 2002, A&A, 385, 337 

Truelove, J. K., Klein, R. I., McKee, C. F. et al. 1997, ApJ, 489, L179 



